蝴蝶优化算法的基本原理¶
蝴蝶优化算法 (Butterfly Optimization Algorithm, BOA）是由印度学者 Sankalap Arora 等人于 2019 年提出的一种新型智能优化算法。该算法通过模拟蝴蝶的受食行为来对最优问题 求解。该算法具有收玫速度快, 寻优能力强的特点。
经过了大自然百万年的自然选择, 蝴蝶凭借种群之间的合作, 以及自身的受食特性得以存活。蝴蝶通过触觉、听觉、嗅觉和视觉来寻找食物、躲避天敌及寻找配偶。在这些感觉中, 嗅觉是最重要的, 蝴蝶的嗅觉感受器几乎覆盖蝴蝶的整个身体, 如触角、腿、触须等。这些嗅觉感受器实际上是蝴蝶体表的神经细胞, 通常也被称为化学感受器。蝴蝶能够借助嗅觉感受器精确地定位某种气味源, 并且每只蝴蝶自身也可以散发出具有一定强度的气味。所以,在蝴蝶优化算法中, 一只蝴蝶对应解空间的一个搜索单位。每只蝴蝶可以释放出一定强度的香味, 这个香味与适应度函数有关, 并且当蝴蝶在解空间内移动时, 其对应的适应度值也会随之变化。在这个过程中, 香味会传播一定的距离, 并且随着蝴蝶之间的距离递减, 其他的 蝴蝶可以感知到这种香味, 这也就是蝴蝶优化算法中的蝴蝶在整个解空间中和其他蝴蝶交换信息的方式。当一只蝴蝶感知到其他蝴蝶身上散发出的香味时, 它会直接向这个散发出香味的蝴蝶移动，这个过程在蝴蝶优化算法中被称为全局搜索; 当蝴蝶感受不到其他蝴蝶身上散发出的香味时, 它会在自己附近做随机探索, 这个过程在蝴蝶优化算法中被称为局部搜索。
蝴蝶的香味¶
在蝴蝶优化算法中, 每只作为搜索单元的蝴蝶都会释放一定强度的香味。每只蝴蝶身上都有自己独特的香味, 并且这些香味的强度会在传播的过程中随着距离递增而逐步衰减, 这一特性也是蝴蝶优化算法和其他启发式算法的主要区别。在蝴蝶优化算法中, 计算香味值 $f$ 与三个重要的参数有关：感知形态 $c$ 、刺激因子 $I$ 及功率指数 $a$ 。
- 感知形态 $c$ 是指感知香味的方式, 在蝴蝶优化算法中, 指代嗅觉。感知形态是算法初始化过程的一个常量, 这也是基础蝴蝶优化算法可优化的一个参数。
- 在蝴蝶优化算法中, 刺激因子 $I$ 是通过当前场景下蝴蝶的适应度函数计算得出的, 所以 当一只蝴蝶的适应度值较大时, 周围的蝴蝶在全局搜索这一过程中就会向适应度值较大的蝴蝶移动。
- 功率指数 $a$ 在蝴蝶优化算法中是一个常数, 它的取值会产生三种效果：线性响应、响应 压缩和响应膨胀。响应膨胀是指随着 $I$ 增大, $f$ 比 $I$ 增长得快, 即 $f$ 对 $I$ 的导函数是增函数; 响应压缩是指随着 $I$ 增大, $f$ 比 $I$ 增长得慢, 即 $f$ 对 $I$ 的导函数是减函数; 线性响应是指 $I$ 和 $f$ 的增长速率相同。通过模拟大自然, 有研究表明, 生物会因环境改变而受到刺激, 但是受到刺激的强度并不是特别强烈。在蝴蝶优化算法中, $a$ 的取值范围为 $(0,1)$ 。
对于蝴蝶优化算法中香味公式 $f$ 的构造, 主要涉及两个问题即 $I$ 的变化和 $f$ 的公式结构。 出于简化模型的考虑, 刺激因子 $I$ 与编码后的目标函数有关。 $f$ 的大小是相对的, 为了从其他感知形态中区分出香味, 需要把 $c$ 作为 $f$ 公式中的一个乘积因子。通过上述对功率指数 $a$ 的分析, 随着 $I$ 的变化, $f$ 的变化不如 $I$ 的变化程度大, 所以 $a$ 应该作为 $I$ 的一个指数因子。通 过上述的讨论, $f$ 为 $$ f=c I^a $$ 其中, $f$ 是每只蝴蝶散发出的香味的衡量值; $c$ 是感知形态; $I$ 是刺激因子; $a$ 是功率指数。 在一般情况下, $a$ 和 $c$ 的取值范围都是 $(0,1)$ 。但在极端的情况下，若 $a=1$, 则表示香味 在传播的过程中没有任何损耗, 也就是说，一只蝴蝶身上散发出的全部香味都会被另一只蝴蝶感受到。在这种理想情况下, 任意一只蝴蝶散发出的香味都会被解空间中的任意一只蝴蝶感受到，这样算法很容易陷入局部最优。若 $a=0$, 则表示一只蝴蝶身上散发出的香味不能被其他任何一只蝴蝶感受到, 也就是香味在传播的路径上全部被消耗掉, 所以 $a$ 的取值对于算法的性能至关重要。另外一个参数 $c$ 影响蝴蝶优化算法的收玫速度。理论上 $c$ 的取值可以是 $(0, \infty)$, 但是实际的取值要与优化问题相结合, 在大多数情况下, $c$ 的取值为 $(0,1)$ 。所以 $a$ 和 $c$ 作为两个常量参数, 决定了整个蝴蝶优化算法的搜索能力和收玫速度。
蝴蝶的移动与迭代¶
关于蝴蝶优化算法的迭代过程, 首先需要约定如下假设:
(1)所有蝴蝶都会散发出一定强度的香味。
(2) 所有蝴蝶要么在原地附近随机搜索, 要么直接向香味值最大的蝴蝶移动。
(3) 刺激因子仅由目标函数决定，不受其他因素影响。
蝴蝶优化算法的执行过程主要分为三个阶段：初始阶段、迭代阶段和终止阶段。每次执行蝴蝶优化算法时, 首先执行初始阶段, 然后根据初始化后的蝴蝶进行迭代, 最后满足终止条件后, 输出结果。
在初始阶段, 首先定义目标函数和解空间, 以及蝴蝶优化算法中的常量, 包括感知形态 $c$ 、功率指数 $a$, 以及切换概率 $p$ 。其次确定蝴蝶的数量, 并随机在解空间中分散所有蝴蝶, 通过蝴蝶优化算法初始化每只蝴蝶, 计算其适应度值和散发出的香味值。
初始阶段后是迭代阶段。在每次迭代中, 蝴蝶在解空间的位置都会重新分布, 所以每只蝴蝶的适应度值和香味值都需要重新计算。在迭代阶段有两个重要的过程: 全局搜索和局部 搜索。在全局搜索中, 蝴蝶朝着香味值最大的蝴蝶 $g^*$ 移动, 即 $$ x_i^{t+1}=x_i^t+\left(r^2 \times g^*-x_i^t\right) \times f_i $$ 其中, $x_i^t$ 表示第 $i$ 只蝴蝶在第 $t$ 次迭代中所对应的解; $g^*$ 表示当前迭代次数中的最优解; $f_i$ 表 示第 $i$ 只蝴蝶散发出的香味值; $r$ 是区间 $[0,1]$ 内的一个随机数。 蝴蚞优化算法的局部搜索公式为 $$ x_i^{t+1}=x_i^t+\left(r^2 \times x_j^i-x_k^t\right) \times f_i $$ 其中, $x_j^t$ 和 $x_k^{\prime}$ 分别表示在第 $t$ 次迭代中第 $j$ 只和第 $k$只蝴蝶所对应的解; $r$是区间 $[0,1]$ 内的随 机数。若 $x_j^{\prime}$ 和 $x_k^{\prime}$ 对应的蝴蝶属于同一个种群, 则式就表示局部的随机移动。
在大自然中, 因为蝴蝶搜索食物包括全局搜索和局部搜索, 所以蝴蝶优化算法也模拟了这两个搜索过程。因为有很多不可预知的天气因素, 包括下雨、大风等, 所以在蝴蝶优化算 法中常用一个常量——切换概率 $p$ 来判断蝴蝶是在做全局搜索还是在做局部搜索, 蝴蝶种群每迭代一次, 每只蝴蝶都会产生一个随机数 $r$, 若当前的随机数 $r>p$, 则进行全局搜索; 否则进行局部搜索。
蝴蝶优化算法的步骤如下:
步煛 1：设置蝴蝶优化算法的基本参数。
步聚 2：随机初始化蝴蝶种群的位置。
步囦 3: 通过适应度函数计算蝴蝶的刺激因子 $I$ 。
步㟫 4: 计算蝴蝶的香味值, 并记录香味值最大的全局最优蝴蝶。
步聚 5：生成随机数 $r$ 并与切换概率 $p$ 进行比较。若 $r>p$, 则利用式进行搜索；否则利用式进行局部搜索。
代码¶
import numpy as np
import random
import numpy as np
from matplotlib import pyplot as plt
import copy
def initialization(pop,ub,lb,dim):
''' 种群初始化函数'''
'''
pop:为种群数量
dim:每个个体的维度
ub:每个维度的变量上边界，维度为[dim,1]
lb:为每个维度的变量下边界，维度为[dim,1]
X:为输出的种群，维度[pop,dim]
'''
X = np.zeros([pop,dim]) #声明空间
for i in range(pop):
for j in range(dim):
X[i,j]=(ub[j]-lb[j])*np.random.random()+lb[j] #生成[lb,ub]之间的随机数
return X
def BorderCheck(X,ub,lb,pop,dim):
'''边界检查函数'''
'''
dim:为每个个体数据的维度大小
X:为输入数据，维度为[pop,dim]
ub:为个体数据上边界，维度为[dim,1]
lb:为个体数据下边界，维度为[dim,1]
pop:为种群数量
'''
for i in range(pop):
for j in range(dim):
if X[i,j]>ub[j]:
X[i,j] = ub[j]
elif X[i,j]<lb[j]:
X[i,j] = lb[j]
return X
def CaculateFitness(X,fun):
'''计算种群的所有个体的适应度值'''
pop = X.shape[0]
fitness = np.zeros([pop, 1])
for i in range(pop):
fitness[i] = fun(X[i, :])
return fitness
def SortFitness(Fit):
'''适应度值排序'''
'''
输入为适应度值
输出为排序后的适应度值，和索引
'''
fitness = np.sort(Fit, axis=0)
index = np.argsort(Fit, axis=0)
return fitness,index
def SortPosition(X,index):
'''根据适应度值对位置进行排序'''
Xnew = np.zeros(X.shape)
for i in range(X.shape[0]):
Xnew[i,:] = X[index[i],:]
return Xnew
def BOA(pop, dim, lb, ub, MaxIter, fun):
'''蝴蝶优化算法'''
'''
输入：
pop:为种群数量
dim:每个个体的维度
ub:为个体上边界信息，维度为[1,dim]
lb:为个体下边界信息，维度为[1,dim]
fun:为适应度函数接口
MaxIter:为最大迭代次数
输出：
GbestScore:最优解对应的适应度值
GbestPositon:最优解
Curve:迭代曲线
'''
p=0.8 #切换概率
power_exponent=0.1 #功率指数a
sensory_modality=0.1 #感知形态c
X = initialization(pop,ub,lb,dim) # 初始化种群
fitness = CaculateFitness(X, fun) # 计算适应度值
indexBest = np.argmin(fitness) #寻找最优适应度位置
GbestScore = fitness[indexBest] #记录最优适应度值
GbestPositon = np.zeros([1,dim])
GbestPositon[0,:] = X[indexBest, :]
X_new = copy.copy(X)
Curve = np.zeros([MaxIter, 1])
for t in range(MaxIter):
print("第"+str(t)+"次迭代")
for i in range(pop):
FP = sensory_modality*(fitness[i]**power_exponent) #刺激强度I的计算
if random.random()<p: #全局搜索
dis = random.random()*random.random()*GbestPositon - X[i,:]
Temp = np.matrix(dis*FP)
X_new[i,:] = X[i,:] + Temp[0,:]
else:#局部搜索
Temp = range(pop)
JK = random.sample(Temp,pop) #随机选择个体
dis=random.random()*random.random()*X[JK[0],:]-X[JK[1],:]
Temp = np.matrix(dis*FP)
X_new[i,:] = X[i,:] + Temp[0,:]
for j in range(dim):
if X_new[i,j] > ub[j]:
X_new[i, j] = ub[j]
if X_new[i,j] < lb[j]:
X_new[i, j] = lb[j]
#如果更优才更新
if(fun(X_new[i,:])<fitness[i]):
X[i,:] = copy.copy(X_new[i,:])
fitness[i] = copy.copy(fun(X_new[i,:]))
X = BorderCheck(X, ub, lb, pop, dim) # 边界检测
fitness = CaculateFitness(X, fun) # 计算适应度值
indexBest = np.argmin(fitness)
if fitness[indexBest] <= GbestScore: # 更新全局最优
GbestScore = copy.copy(fitness[indexBest])
GbestPositon[0,:] = copy.copy(X[indexBest, :])
Curve[t] = GbestScore
return GbestScore, GbestPositon, Curve
'''适应度函数'''
def fun(X):
O=X[0]**2 + X[1]**2
return O
'''蝴蝶优化算法求解x1^2 + x2^2的最小值'''
'''主函数 '''
#设置参数
pop = 50 #种群数量
MaxIter = 100 #最大迭代次数
dim = 2 #维度
lb = -10*np.ones(dim) #下边界
ub = 10*np.ones(dim)#上边界
#适应度函数选择
fobj = fun
GbestScore,GbestPositon,Curve = BOA(pop,dim,lb,ub,MaxIter,fobj)
print('最优适应度值：',GbestScore)
print('最优解[x1,x2]：',GbestPositon)
#绘制适应度曲线
plt.figure(1)
plt.plot(Curve,'r-',linewidth=2)
plt.xlabel('Iteration',fontsize='medium')
plt.ylabel("Fitness",fontsize='medium')
plt.grid()
plt.title('BOA',fontsize='large')
plt.show()
参考资料：
- 范旭 编 《Python智能优化算法：从原理到代码实现与应用》
